#include<fstream>
#include<iostream>
#include<complex>
#include<vector>
#include<cmath>
#include<bitset>
#include"system_ho.h"
//#include<boost/math/special_functions/bessel.hpp>
#include"help.h"
using namespace std;
//using namespace boost::math;

typedef complex<double> complexd;

// double Vr(double r)
// {
//   //  return 1.;
//   return exp(-r)/r;
// }

// const int nmax=20;
// const int points=100;
// const double rmax=20;
// const double lengthb=1.;
// const int l=0;
// vector<double> x(points),w(points);
// double Vrho(int a,int b)
// {
//   double res(0);
//   for(int i=0;i<points;i++)
//     {
//       res+=ho_r(a,l,lengthb,x[i])*Vr(x[i])*ho_r(b,l,lengthb,x[i])*x[i]*x[i]*w[i];
//     }
//   return res;
// }
// double newVrr(double r1,double r2)
// {
//   double res(0);
//   for(int a=0;a<nmax;a++)
//     for(int b=0;b<nmax;b++)
//       {
// 	res+=ho_r(a,l,lengthb,r1)*ho_r(b,l,lengthb,r2)*Vrho(a,b);
//       }
//   return res;
// }
// double newVr(double r)
// {
//   double res(0);
//   for(int i=0;i<points;i++)
//     {
//       res+=newVrr(r,x[i])*x[i]*x[i]*w[i];
//     }
//   return res;
// }
// template<size_t N>
// bitset<N> operator + (bitset<N> a, bitset<N> b)
// {
//   bitset<N> ans;
//   while(b.any())
//     {
//       ans = a^b;
//       b = (a&b)<<1;
//       a = ans;
//     }
//   return ans;
// }
// template<size_t N>
// bitset<N> operator - (const bitset<N>&a)
// {
//   return ~a + bitset<N>(1);
// }
// template<size_t N>
// bitset<N> operator - (bitset<N> a,bitset<N> b)
// {
//   return a + (-b);
// }

// template<size_t N>
// int ctz(bitset<N> a)
// {
//   if(a.none()) return 0;
//   int count=0;
//   while(!a[0])
//     {
//       a>>=1;
//       ++count;
//     }
//   return count;
// }
// template<size_t N>
// bitset<N> next(const bitset<N>&a)
// {
//   bitset<N> one(1);
//   bitset<N> t= a | (a-one);
//   return (t + one) | (((~t & -~t) - one) >> (ctz(a) + 1));
// }


int main()
{
  // gauleg(0,rmax,x,w);
  // for(int i=1;i<100;i++)
  //   {
  //     double r=0.1*i;
  //     cout<<r<<"\t"<<Vr(r)<<"\t"<<newVr(r)<<endl;
  //   }
  // int a;
  // cin>>a;
  // bitset<8> b(a);
  // cout<<b<<endl;
  // cout<<next(b)<<endl;
  // b<<=2;
  // bitset<8> one(1);
  // b=b-one;
  // cout<<b<<endl;
  // const int nmax=100;
  // const double step=0.1;
  // const double theta=15;
  // for(int nx=0;nx<nmax;nx++)
  //   //    for(int ny=0;ny<nmax;ny++)
  //     {
	
  // 	complexd v=cmplx_ho_k(0,1,1.,nx*step*complexd(cos(theta),sin(theta)));
  // 	//	cout<<nx*step<<"\t"<<ny*step<<"\t"<<v.real()<<"\t"<<v.imag()<<endl;
  // 	cout<<nx*step<<"\t"<<v.real()<<"\t"<<v.imag()<<endl;	
  //     }
  System_Tz_HO ho("./sp.dat","./vint.dat");
  ho.setup();
  ho.print2B();
  return 0;
}
